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1.  Summary 


This  work  reports  anisotropic  and  isotropic  elastic  properties  of  dislocations  and  cracks  in  the 
a-phase  of  the  energetic  molecular  crystal  hcxahydro- 1 ,3,5-trinitro-s-triazine  (aRDX).  The 
properties  determined  in  this  work  can  be  used  to  determine  probable  slip  systems  in  aRDX  and 
the  structure  of  dislocations.  The  procedures  given  in  this  work  for  determination  of  the 
anisotropic  dislocation  line  energy  and  crack  tip  dislocation  nucleation  load  factors  follow  those 
outlined  by  Knap  ( 1 ).  This  work  contains  a  summary  of  the  anisotropic  elastic  solutions  for 
infinite  straight  dislocations  given  by  Barnett  and  Swanger  (2)  and  Barnett  and  Asaro  (3)  and  the 
isotropic  and  anisotropic  dislocation  nucleation  models  of  Rice  (4)  and  Sun  and  Beltz  (5).  The 
appendix  contains  the  Matlab  implementation  of  the  algorithms  given  in  the  text. 


2.  Elastic  Constants 


The  orthotropic  elastic  constants  used  for  aRDX  are  given  in  this  section.  The  elastic  constants 
are  given  in  the  material  coordinates,  X„  which  are  aligned  with  the  major  axes  of  the  material 
stiffness  tensor.  The  linear  elastic  relationship  between  the  stress  and  strain  tensors,  a  and  e,  are 
related  by  a  fourth  order  tensor,  Cp/.  Voigt  form  uses  major  and  minor  symmetry  of  the  stress 
and  strain  tensors  to  write  the  fourth-order  tensor,  Cp/,  as  a  6  x  6  matrix,  cy.  For  an  orthotropic 
material,  the  stress-strain  relationship  in  Voigt  form  is  given  by 
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Previously  reported  (<5,  7)  molecular  dynamics  simulations  under  an  applied  strain  were  used  to 
determine  the  components  of  cy.  Each  volumetric  strain  component,  results  in  a  triaxial  state 
of  stress,  a#,  due  to  Poisson’s  effect  and  is  used  to  determine  the  volumetric  components  of  cy. 
Each  shear  strain  results  in  a  single  linearly  proportional  shear  stress,  giving  the  shear 
components  of  cy.  The  compliance  tensor,  S p/,  is  given  by  the  inverse  of  the  stiffness  tensor, 

Cpz,  but  there  is  not  a  method  available  for  taking  the  inverse  of  a  fourth-order  tensor  and,  due  to 
material  symmetries,  the  81  components  of  Cp/  would  produce  a  rank  deficient  9x9  matrix.  The 
Voigt  form  of  the  stiffness  tensor,  cy,  includes  the  crystal  symmetries  and  therefore  is  not  rank 
deficient  and  an  inverse  exists.  The  Voigt  form  of  the  compliance  tensor,  s, y,  given  by  the  inverse 
of  cy,  in  terms  of  common  orthotropic  engineering  terms  is  given  by 
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where  E,  is  the  orthotropic  Young’s  modulus,  vy  is  the  orthotropic  Poisson’s  ratio  and  Gy  is  the 
orthotropic  shear  modulus.  Values  for  cy,  sy,  Eh  vy,  and  Gy  are  given  in  table  1. 

Table  1.  Orthotropic  elastic  and  engineering  constants 
for  aRDX  (6,  7). 


Stiffness 

Components 

(GPa) 

Compliance 

Components 

(GPa1) 

Engineering 
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(GPa) 
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Voigt  and  Reuss  averages  are  used  to  find  the  isotropic  bulk  and  shear  modulii.  The  Voigt 
average  assumes  a  uniform  strain  field  leading  to  an  overestimate  of  the  stress  field  calculated 
from  anistropic  elasticity.  The  Reuss  average  assumes  uniform  stress  and  results  in  an 
underestimate  of  the  stress  field  from  anisotropic  elasticity.  The  Voigt  average  for  the  bulk 
modulus,  Bv ,  and  shear  modulus,  Gv,  found  from  the  stiffness  tensor  are  given  by 


Bv  = 


i,j 


-( r „ 


(3) 


and  are  tabulated  in  table  2  for  aRDX.  Likewise,  the  Reuss  average  for  the  bulk  modulus,  Br, 
and  shear  modulus,  Gr,  found  from  the  compliance  tensor  are  given  by 
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and  are  also  tabulated  in  table  2  for  aRDX.  Table  2  also  presents  isotropic  averages  for  the 
Young’s  modulus  and  Poisson  ratio  found  for  either  Voigt  or  Reuss  averages  of  the  bulk  and 
shear  modulus  using  the  following  relationships 
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Table  2.  Isotropic  averages  of 

orthotropic  elastic  constants 
for  aRDX. 


Voigt 

(GPa) 

Reuss 

(GPa) 

Bv 

14.0 

Br 

14.0 

Gv 

6.2 

Gr 

5.5 

Ev 

16.2 

Er 

14.5 

vv 

0.31 

Vr 

0.33 

The  elastic  constants  presented  in  table  1  are  given  in  the  laboratory  or  material  coordinates,  Xt, 
which  are  the  coordinates  where  the  stiffness  and  compliance  tensors  from  equations  1  and  2  are 
in  their  simplest  forms,  showing  the  symmetry  of  the  crystal.  The  material  coordinates,  Xj,  are 
normally  not  coincident  with  the  crack  coordinates  x,  as  shown  in  figure  1.  The  crack  coordinates 
shown  in  figure  1  are  aligned  so  that  the  X2  direction  is  normal  to  the  crack  face  and  the  x/ 
direction  is  aligned  with  the  crack  front.  The  dislocation  geometry  is  also  described  in  terms  of 
these  coordinates  where  the  dislocation  line  direction  of  the  dislocation  to  be  emitted  from  the 
crack  tip  will  be  in  the  X3  direction. 

The  anisotropic  elasticity  equations  in  the  following  sections  were  derived  by  Barnett  and 
Swanger  (2)  to  yield  expressions  that  utilize  the  material  coordinates.  The  equations  also  utilize 
the  fourth  order  tensor  version  of  the  stiffness  tensor,  Cijkl.  An  algorithm  for  creating  a  fourth 
order  tensor  from  the  6x6  Voigt  form  is  given  in  the  appendix. 
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Figure  1.  Coordinates  used  in  elastic  crack  tip  solution.  X,  are  the  material 

coordinates  aligning  with  the  major  axis  of  the  material  stiffness  tensor, 
Xj  are  the  crack  tip  coordinates,  where  xt  is  the  direction  of  the  crack 
front,  x2  is  normal  to  the  crack  face,  and  x3  is  the  plane  dimension.  K,  are 
the  applied  load  intensity  factors  describing  the  local  modes  of 
deformation  at  the  crack  tip. 


3.  Anisotropic  Elastic  Dislocation  Line  Energy 


This  section  will  give  a  brief  description  of  how  to  determine  the  dislocation  line  energy  using 
the  methodology  developed  by  Barnett  and  Swanger  (2).  This  section  is  a  reproduction  of  their 
work  and  will  therefore  use  the  same  constants  and  terminology.  In  their  work  they  refer  to  the 
coordinates  in  which  the  elastic  constants  are  defined  in  as  the  laboratory  or  material  coordinates, 
Xj.  In  equation  1  of  Barnett  and  Swanger  (2),  the  energy  per  unit  length  of  a  dislocation  line  is 
given  by 

R 

E  =  Kmgbmbg\n—  where  (jn,g  =  1,2,3)  (6) 

r0 

where  Kmg=Kgm  is  the  energy  prefactor  that  only  depends  on  the  elastic  constants  C,^/  and  the 
orientation  of  the  dislocation  with  respect  to  the  laboratory  Xj  coordinates.  The  other  parameters 
in  equation  6  are  as  follows  b  is  the  burgers  vector,  ra  is  the  radius  of  the  dislocation  core,  R  is  a 
radial  dimension  on  the  same  order  as  the  specimen.  In  equation  AI.l  of  appendix  I  of  Barnett 
and  Asaro  (3),  a  straightforward  method  for  evaluating  Kmg  is  given  as 


K mg  Qn2  £PJW 
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ngip^wmrs 
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where  is  the  alternating  or  Levi-Cevita  tensor,  t  is  the  dislocation  line  direction,  z  is  any  unit 
vector  perpendicular  to  t,  and  t//  is  an  angular  variable  in  the  plane  z*t=0  (see  figure  2).  Mir  is  the 
symmetric  Christoffel  stiffness  matrix  given  by  Mir  =  CjjrsZjZs ■  M[rl  is  the  inverse  of  Mjr. 
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X*  [001] 


Figure  2.  Spherical  polar  coordinates  used  to  define  the 
dislocation  line  direction,  t,  with  respect  to  the 
crystal  axes  or  laboratory  axes,  X z  is  the  unit 
normal  vector  to  t  and  \|/  is  the  angle  of  rotation 
integrated  over  in  equation  7. 


In  our  work,  the  dislocation  line  direction,  t,  with  respect  to  the  laboratory  axis,  Xj,  is  known  a 
priori  and  is  fixed.  Using  the  spherical  polar  coordinates  shown  in  figure  2,  t  is  given  by 

(sirup  cos  St 

sin  cp  sin  0  |  (8) 

cos  cp  ) 

With  t  fixed,  any  unit  vector  z  normal  to  t  and  its  derivative  are  given  by 
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The  integral  in  equation  7  is  evaluated  using  trapezoidal  numerical  integration  yielding  a  fourth 
order  tensor  that  is  evaluated  with  the  remaining  terms  of  equation  7  yielding  Kmg. 

In  table  3,  values  for  the  anisotropic  prelogarithmic  factor,  Kmgbmbg.  are  given  for  several  slip 
systems  reported  on  by  Munday  (7),  Munday  et  al.  (5),  and  Mathew  et  al.  (9).  Isotropic  values 
are  also  given  using  the  Voigt  and  Reuss  averages.  The  dislocation  line  energy  prefactor  for  an 
isotropic  material  is  given  by 
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(11) 


Kiso 


G(1  —  v  cos2  6 ) 
4n(l  —  v) 


where  G  and  v  are  the  Voight  or  Reuss  isotropic  averages  for  the  shear  modulus  and  Poisson’s 
ratio,  respectively,  and  9  is  the  angle  between  the  line  direction,  t,  and  Burger’s  vector,  b.  Plots 
of  the  prelogarithmic  energy  factors  as  a  function  9  are  shown  in  figure  3  for  b  =  Vi[  100] 
dislocations  on  the  (010)  or  (001)  glide  planes.  Matlab  algorithms  for  determining  the  line 
energies  are  given  in  the  appendix. 

Table  3.  Anisotropic  and  isotropic  prelogarithmic  energy  factor  ( E  =  Kmgbmbg)  for  aRDX  slip 
systems. 


Slip  plane, 

n 

Burgers 
vector,  b 

Line 

direction,  t 

Dislocation 

type 

E  Aniso 
(J  nrhe-O 

E  Voigt 
(J  m 1)e-9 

E  Reuss 
(J  m 1)e-9 

(010) 
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[001] 
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0.32 
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screw 
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0.22 

0.19 

(010) 

[001] 

[100] 

edge 
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screw 
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0.47 

(001) 
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screw 
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0.56 

(001) 
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screw 

0.22 

0.22 
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(001) 

y2[ioo] 

[010] 

edge 

0.33 
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y2[oio] 

[100] 

edge 

0.19 

0.23 

0.21 

(021) 

y>[ol2] 

[012] 

screw 

0.64 

0.69 

0.61 

(021) 

[012] 

[012] 

screw 

2.56 

2.75 

2.42 

(021) 

y2[ioo] 

[012] 

edge 

0.36 

0.32 

0.29 

(021) 

y2[ioo] 

[100] 

screw 

0.22 

0.22 

0.19 

(021) 

[100] 

[100] 

screw 

0.9 

0.88 

0.78 

(011) 

[100] 

[100] 

screw 

0.9 

0.88 

0.78 

(011) 

[100] 

[100] 

edge 

1.37 

1.27 

1.16 

(011) 

[011] 

[011] 

screw 

1.28 

1.17 

1.03 

'  ■  ■'  !  '  Mi .  I  "  ! :  jj  uiv 


Figure  3.  Anisotropic  and  isotropic  energy  prelog  factors  as  a  function  of  dislocation  type  for  dislocations 
with  b  =  Vi\\  0  0]  on  the  (a)  (010)  and  (b)  (001)  glide  planes. 
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4.  Anisotropic  Dislocation  Nucleation 


A  summary  of  the  anisotropic  version  of  Rice’s  ( 4 )  dislocation  nucleation  continuum  model 
developed  by  Sun  and  Beltz  (5)  will  be  presented  followed  by  simplifications  made  for  the 
systems  under  pure  mode  II  or  mode  III  loading.  In  this  work,  the  dislocation  line  direction,  t, 
given  in  equation  8  and  shown  in  figure  2  coincides  with  the  crack  front,  given  as  the 
Vi-coordinate  in  figure  4. 


Figure  4.  Crack  tip  geometry  where  the  red  plane  is  coincident  with  the  crack 
plane,  the  crack  front  coincides  with  the  .^-coordinate  and  the 
dislocation  line  direction,  t,  and  the  slip  plane  is  shown  in  blue  and 
is  tilted  by  an  angle  0  relative  to  the  crack  plane.  The  slip  direction 
is  given  by  5  and  angle  4>. 


In  equation  58  of  Barnett  and  Asaro  (3),  the  Griffith  energy  for  a  loaded  crack  tip  in  an 
anisotropic  material  is 

1  _± 

G  —  ^kmKmgkg 


(12) 


where  Kmg  is  the  anisotropic  energy  prefactor  given  by  equation  7,  and  k,„  =  (K1,K2,K3) 

=  (KII,KI,KIII)  are  the  applied  stress  intensity  factors  as  shown  in  figure  4.  In  fracture 
mechanics,  equation  12  is  normally  given  as 

G  ~  km^mgkg  (13) 

where  Amg  =  K^g/8n.  Kmg  and  Amg  are  symmetric  and  can  have  off-diagonal  elements. 


Sun  and  Beltz  (5)  give  the  modified  dislocation  nucleation  criterion  of  Rice  ( 4 )  for  an  anisotropic 
media  as 


(14) 
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where  sm(<j>)  =  [cos  (/>.  0,  sin  (f>]  is  the  orientation  of  the  slip  direction  5  on  the  slip  plane  and  yus  is 

6  f  f 

the  unstable  stacking  fault  energy  found  from  the  generalized  stacking  fault  energy  surface.  /cm 
is  the  effective  stress  intensity  factor  on  the  slip  plane  tilted  at  an  angle  of  9  relative  to  the  crack 
plane  and  is  given  by  Rice  ( 4 )  as 


kmf  =  Pmgi.S)kg  = 


cos- 

2 


3  sin2  - 
2 


2  9  .  9  n 

cos  -sin-  0 

2  2 


o  2  6  ■  t) 

-3cosz  -sin- 
2  2 


e 


3  0 

cos J  - 
2 

0 


0 


cos-j 

2J 


(15) 


Sun  and  Beltz  (5)  point  out  that  the  term  smk ^  is  the  effective  stress  intensity  factor  due  to  the 
combination  of  KII  and  Kill  loading  along  the  sm  direction  of  the  slip  plane. 


The  function  p((|),0)  describes  the  effective  material  properties  in  the  direction  of  slip  and  is  given 
by 


P(0>  9)  =  8nsmK^sg 


(16) 


where  8nK„,g  from  equation  12  is  used  instead  of  Amg.  Kmg  is  related  to  Kmg  given  in  equation  7 
by  the  rotation  from  the  laboratory  coordinates,  Xt,  to  the  crack  coordinates,  in  figure  4,  given 
by 


A mg 


RmiKijRjg 


(17) 


The  rotation  matrix,  Rg,  rotating  Km„  to  the  crack  coordinates,  Xj,  by  0  about  the  xj  axis  is 


cos  0 
—  sin  9 
.  0 


sin  6 
cos  6 
0 


O' 

0 

1. 


(18) 


In  this  work,  the  crack  tip  geometry  is  simplified  as  shown  in  figure  5,  where  the  slip  plane  is 
coincident  with  the  crack  plane,  9  =  0.  The  rotation  matrix  given  by  equation  18  is  equal  to  the 
identity  and  equation  17  simplifies  to  K-  -  1  =  Kij.  For  the  emission  of  a  pure  edge  dislocation, 
the  slip  direction,  5,  coincides  with  the  xj  or  KII  direction,  where  (|)  =  0  and  5„,(0)  =  [100],  and 
equation  16  simplifies  to  p(0,0)  =  StiA'i  i  .  The  stress  intensity  factors  for  this  configuration  given 
by  equation  14  are 

KII  =  j8nKllYus  (19) 

where  the  Kmg  tensor  with  component  K\  \  is  found  from  equation  7  for  t  =  [KII  0  0]  and  y„v  is  the 
energy  barrier  to  slip  in  the  5  direction. 
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Figure  5.  Simplified  crack  tip  geometry  for  the  emission  of  a  pure  edge 
dislocation  where  the  blue  slip  plane  is  coincident  with  the  red 
crack  plane  (0  =  0)  and  the  slip  direction,  5,  is  aligned  with  the  Kll 
direction  (<|>  =  0). 

Alternatively,  for  the  emission  of  a  pure  screw  dislocation  shown  in  figure  6,  the  slip  direction,  5, 
coincides  with  the  xj  or  Kill  direction  where  (j)  =  jt/2,  sm( n/2)  =  [001],  and  equation  16  simplifies 
to  p{ jt/2,0)  =  %nKy]-  The  stress  intensity  factors  for  this  configuration  is  given  by  equation  14  is 

Kill  =  j8nK33Yus  (20) 

where  K 33  component  now  to  a  Kmg  tensor  given  by  equation  7  for  t  =  [0  0  Kill]  and  yus  is  the 
energy  barrier  to  slip  in  the  5  direction.  Stress  intensity  factors  for  isotropic  averaged  material 
properties  are  given  by  substituting  in  equations  19  or  20  values  for  K  determined  from 
equation  11. 


Figure  6.  Simplified  crack  tip  geometry  for  the  emission  of  a  pure  screw 
dislocation  where  the  blue  slip  plane  is  coincident  with  the  red 
crack  plane  (0  =  0)  and  the  slip  direction,  5,  is  aligned  with  the  Kill 
direction  (c[>  =  Jt/2). 
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Generalized  stacking  fault  (GSF)  energy  surfaces  for  two  aRDX  slip  systems  are  shown  in 
figure  7  for  the  (010)  slip  plane  (7,  8).  These  GSF  surfaces  provide  energy  barriers,  yus,  for  a 
crack  tip  geometry  with  an  X2  =  [010]  crack  face  and  xj  =  t=[001]  or  [100]  crack  front 
/dislocation  hne  direction.  Anisotropic  and  isotropic  critical  stress  intensity  factors  for  these 
crack  tip  geometries  from  equations  19  and  20  are  given  in  table  4  along  with  values  for  Kmg 
from  equation  7  and  11.  Using  anisotropic  material  properties,  the  lowest  nucleation  threshold 
for  the  (010)  crack  face  is  a  screw  type  ViflOO]  dislocation  emitted  under  Kill  loading.  This  slip 
system  also  had  the  lowest  dislocation  line  energy  in  table  3. 


Figure  7.  GSF  energy  surfaces  for  dislocations  on  the  (010)  slip  plane. 


Table  4.  Anisotropic  and  isotorpic  critical  stress  intensity  factors  for  the  (010)  crack  face  in  aRDX. 


Slip  plane 

8 

/,  line  dir 

Dislocation  Type/ 
Load 

Yus  o 
(mJ/m 

Aniso 

Kmg(UMm) 

Voigt  Reuss 

Critical  AT/ or  Kill  (Mpa  rrA2) 
Aniso  Voigt  Reuss 

(010) 

[100] 

[001] 

edge  /  K1I 

164 

826 

714 

647 

0.058 

0.054 

0.052 

(010) 

[001] 

[001] 

screw  /  Kill 

351 

320 

495 

436 

0.053 

0.045 

0.042 

(010) 

[100] 

[100] 

screw  /  Kill 

164 

576 

495 

436 

0.049 

0.066 

0.062 

(010) 

[001] 

[100] 

edge  /  K1I 

351 

503 

714 

647 

0.067 

0.079 

0.076 

The  (010)[100]  GSF  surface  in  figure  7  has  a  local  energy  minimum  indicating  a  stable  stacking 
fault  configuration  for  this  slip  system.  In  the  case  of  dislocation  nucleation,  a  stable  stacking 
fault  means  a  partial  dislocation  will  be  emitted  from  the  crack  tip.  The  shear  stress  at  the  crack 
front  will  cause  the  partial  dislocation  to  move  away  from  tip,  leaving  behind  a  stacking  fault. 
The  stacking  fault  material  is  at  a  higher  energy  pulling  the  partial  dislocation  back  toward  the 
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crack  tip.  There  is  also  an  image  force  pulling  the  partial  dislocation  back  toward  the  crack  tip.  A 
stable  equilibrium  position,  rA,  for  the  partial  dislocation  is  found  when  the  Peach  Koehler  forces 
due  to  the  externally  loaded  crack  tip  are  balanced  out  by  the  stacking  fault  and  image  forces, 
leading  to  equation  22  in  Sun  and  Beltz  (5)  and  reproduced  here  to  include  lattice  frictional 
forces 


kAbA 

pnrA 


( Ysf  +  0PbA) 


8nrA 


(21) 


where  kA  is  the  applied  stress  intensity  factor,  bA  is  the  Burgers  vector  of  the  partial  dislocation 
being  emitted,  p(<f>A)  is  from  equation  16,  and  apbA  is  the  friction  force  due  to  lattice  resistance 
where  op  is  the  Peierls  stress.  For  the  edge-type  partial  dislocations  described  in  table  4,  ka  =  KII 
and  (f>A  =  0,  leading  to 


KlIbA 
V 2nrA 


(: Ysf  +  °pbA) 


h  fonK]  j 

rA 


=  0 


(22) 


The  equilibrium  position,  ra,  is  found  using  equation  24  in  Sun  and  Beltz,  reproduced  here  for 
the  edge-type  partial  dislocations 


where  KII  is  the  load  being  applied,  KIIAniso  is  the  critical  anisotropic  load  factor  calculated  for 
emission  of  the  partial  dislocation  and  KII  >KlIAniso ,  meaning  the  partial  has  been  nucleated.  The 
isotropic  version  of  equation  23  is  given  by  equation  38  in  Rice  ( 4 )  and  reproduced  here  for  an 
edge  dislocation  with  the  addition  of  lattice  friction 


Y /  cn  — 


(  KII  \2  (Gbj) 

\K 1 1 crit-'  V  Yus  ) 


47t(1  —  v ) 


1  1 

li  ( KU  V 

■  / Ysf  +  <rpbA\ 

\ 

yKiiiS0) 

V  Yus  ) 

(24) 


where  KIIAniso  is  the  isotropic  load  factor  calculated  for  emission  of  the  partial  dislocation. 

Lattice  friction  is  usually  not  included  for  most  materials  because  the  Peierls  stress,  op,  is  usually 
on  the  order  of  10'  to  10  'G.  The  lattice  friction  cannot  be  ignored  for  aRDX  because  the 
Burger’s  vectors  are  large  and  Mathew  et  al.  (9)  found  ap  to  be  approximately  3%-10%  of  G. 
The  large  lattice  friction  could  be  due  to  steric  interactions.  Values  for  rA  and  riso  for  the 
(010)1/2[100]  slip  system  are  given  in  table  5  for  both  with  and  without  lattice  friction.  For  cases 
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with  lattice  friction,  alternative  values  for  op  are  presented.  The  “Iso  Fit”  and  “Aniso  Fit”  op 
values  were  found  by  fitting  to  r  from  the  simulation  data  for  KII  =  0.07  MPa  m  2.  The  MD  and 
GSF  op  values  were  found  by  Mathew  et  al.  (9)  using  molecular  dynamics  and  the  Peierls 
Nabarro  model  with  the  GSF  data  from  Munday  (8)  and  shown  in  figure  7.  Including  op  in 
equations  23  and  24  produces  imaginary  numbers  for  several  values  of  KII.  The  values  produced 
by  the  molecular  dynamics  simulations  (10)  of  a  crack  tip  given  in  the  “Sim.”  column  of  table  5 
do  not  match  any  of  the  op  values  over  the  range  of  KII  loads. 

Table  5.  Location  of  the  first  (OlO)1/^!  100]  partial  edge  dislocation  emitted  from  the  crack  tip  for  the  anisotropic 

and  isotropic  models.  Each  row  corresponds  to  the  ap  value  given  at  the  top  of  the  table.  The  “Sim.”  shows 
values  from  a  molecular  dynamics  simulation  (10). 


op  (GPa) 

ISOTROPIC  MODEL 

ANISOTROPIC  MODEL 

Iso  Fit 

Aniso  Fit 

GSF 

MD 

None 

130 

None 

115 

201 

232 

KII  (Mpa  m'/2) 

Sim.  (A) 

fiso(A) 

Hso(A) 

rA(A) 

rA(A) 

rA(A) 

rA(A) 

0.054 

135 

NaN 

NaN 

NaN 

NaN 

NaN 

0.058 

175 

25 

156 

NaN 

NaN 

NaN 

0.060 

189 

31 

170 

NaN 

NaN 

NaN 

0.065 

235 

47 

215 

44 

NaN 

NaN 

0.070 

62 

284 

63 

263 

62 

17 

NaN 

0.075 

336 

79 

315 

80 

34 

23 

0.080 

392 

95 

369 

98 

45 

35 

0.085 

452 

113 

427 

117 

57 

45 

0.090 

76 

514 

131 

489 

137 

69 

55 

After  the  partial  dislocation  is  nucleated  the  material  at  the  crack  tip  is  in  the  stacking  fault 
configuration.  The  increased  energy  of  the  stacking  fault  configuration  lowers  the  energy  barrier 
for  emission  of  the  second  partial  dislocation.  For  the  case  where  the  partial  dislocation’s 
Burgers  vectors  are  in  the  same  direction  as  the  full  dislocation,  the  first  partial  dislocation 
completely  shields  the  crack  tip  from  the  kcru  load  at  which  it  was  emitted.  For  the  edge 
dislocation,  the  load  for  emission  of  the  second  partial  dislocation  is  then  reduced  to 

KII  =  j87TKn(yus-ys/)  (25) 

Values  for  critical  KII  loads  for  emission  of  the  first  and  second  partial  dislocation  are  given  in 
table  6  for  several  slip  systems  that  were  determined  by  Munday  et  al.  (8)  to  produce  partial 
dislocations  in  aRDX. 
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Table  6.  Anisotropic  and  isotropic  critical  stress  intensity  factors  for  the  (010)  crack  face  in  aRDX. 
Unstable  stacking  fault  energy,  yus,  given  by  Munday  et  al.  (8). 


Slip 

plane, 

n 

Burgers 

vector, 

b 

Line 

direction, 

t 

Type 

Yus 

(mJ  nr2) 

KII  Aniso 

KII  Voigt 

KII  Reuss 

1st 

2nd 

1st 

2nd 

1st 

2nd 

(010) 

wool 

[001] 

edge 

164 

0.058 

0.036 

0.054 

0.034 

0.052 

0.032 

(001) 

wool 

[0101 

edge 

260 

0.070 

0.032 

0.068 

0.031 

0.065 

0.030 

(011) 

wool 

[011] 

edge 

255 

0.070 

0.047 

0.068 

0.045 

0.064 

0.043 

(021) 

wool 

[012] 

edge 

250 

0.071 

0.036 

0.067 

0.034 

0.064 

0.032 

After  the  second  partial  dislocation  is  emitted,  the  two  partial  dislocations  remain  separated  by 
the  stacking  fault.  Even  in  the  absence  of  a  defect,  it  is  energetically  favorable  for  a  full 
dislocation  to  dissosciate  into  two  partials  with  a  fractional  Burgers  vector  separated  by  a 
stacking  fault  because  the  dislocation  line  energy  given  by  equation  6  is  proportional  to  l/?l“.  In 
the  bulk  crystal  the  partial  dislocations  are  pushed  apart  by  elastic  forces  due  to  both  partial 
dislocations  having  similar  Burger’s  vectors.  As  the  partial  dislocations  move  away  from  one 
another  they  increase  the  size  of  the  stacking  fault  resulting  in  an  attractive  force  holding  them 
together.  An  equilibrium  configuration  is  reached  where  the  two  forces  are  balanced  and  are 
described  in  equations  10-14  and  13-150  of  Hirth  and  Lothe  (11).  For  the  pure  edge-type  partial 
dislocations  given  in  table  6,  the  equilibrium  equation  from  Hirth  and  Lothe  is  simplified  as 

2Kmgbm.bg 

ra  = - - -  (2< 

Ysf 

Values  for  ra  are  given  in  table  7  for  anisotropic  and  isotropic  Kmg  values. 

Table  7.  Anisotropic  and  isotorpic  partial  dislocation  separation  distances  for  b  =  x/^[  100]  edge-type 
dislocations. 


Slip 

plane,  n 

Burgers 
vector,  b 

Line 

direction, 

t 

Dislocation 

type 

Ysf 

(mJ  nr2) 

ra(A) 

Aniso 

Voigt 

Reuss 

(010) 

W00] 

[001] 

edge 

101 

73 

63 

65 

(001) 

W00] 

[010] 

edge 

206 

32 

31 

32 

(011) 

W00] 

[011] 

edge 

140 

49 

46 

47 

(021) 

W00] 

[012] 

edge 

187 

38 

34 

35 
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5.  Conclusion 


Several  dislocation  and  crack  tip  properties  were  presented  in  this  work  that  are  useful  for 
energetically  ranking  slip  system  activity  and  predicting  the  dislocation  structure.  Both  isotropic 
and  anisotropic  results  were  provided  and  the  isotropic  results  are  shown  to  predict  dislocation 
nucleation  at  smaller  thresholds  than  the  anisotropic  model.  The  isotropic  and  anisotropic  models 
provided  the  same  energetic  ranking  based  on  dislocation  line  direction.  For  isotropic  materials, 
the  line  energy  of  a  screw  dislocation  is  always  lower  than  that  of  edge  dislocations  due  to  edge 
dislocations  being  scaled  by  (1-n)"1  but  for  anisotropic  materials  this  is  not  always  the  case,  see 
for  instance  NiAl  (12)  and  aFe  (13).  However,  for  the  slip  systems  we  studied  in  this  work,  the 
anisotropic  line  energy  for  screw  dislocations  was  always  lower  than  edge  dislocations  and  was 
in  general  agreement  with  the  isotropic  model. 
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Appendix.  Matlab  Code  for  Dislocation  Line  Energies  and  Nucleation 


This  appendix  appears  in  its  original  form,  without  editorial  change. 
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%  Determine  Anistropic  and  isotropic  dislocation  line  energy  and 
%  dislocation  nucleation  form  a  crack  tip. 

"6 

%  This  script  provides  the  data  for  the  Anisotropic  dislocation  line  energy 
%  and  crack  tip  dislocation  nucleation  ARL  technical  report.  Equation 
%  numbers  given  in  this  script  correspond  to  equation  numbers  given  in  the 
%  report. 

"6 

%  Anistropic  dislocation  orientation  prefactor  tensor  determined  from 
%  equations  given  by  Barnett  &  Swanger  ( PhysStat . Sol  1971)  and  Barnett  & 

%  Asaro  (JMPS  1972)  . 

"6 

%  KII  isotropic  and  anisotropic  crack  tip  loading  factors  for  dislocation 
%  nucleation  using  in  Rice's  (JMPS  1994)  isotropic  dislocation  nucleation 
%  model  or  Sun  &  Beltz's  (JMPS  1994)  anisotropic  dislocation  nucleation 
%  model. 

"6 

%  Lynn  Munday,  August  2012 
"6 

close  all  hidden, clear  all,  clc 
format  short 


% - MANUAL  INPUT - 

"6 

%  RDX  stable  and  unstable  stacking  energy  J/mA2=Pa*m 
%  From  Munday  thesis  (UMd  2011)  and  Munday  et  al .  (Phil  Mag  2012) 
y_us=l 64e-3 ;  %  for  (010) [100] 

y_sf=101e-3; 


%  y_us=260e-3; 
%  y_sf=206e-3; 
"6 

%  y_us=255e-3; 
%  y_sf=140e-3; 
"6 

%  y_us=250e-3; 
%  y_sf=187e-3; 


for  (001) [100] 

for  (011) [100] 

for  (021) [100] 


%lattice  vectors  a,b,c  or  [100]  [010]  [001] 

latVec= [13 . 366  0  0;... 

0  11 . 334  0;  .  .  . 

0  0  10.341] *le-10; 


bvec_n=0 . 5* [ 1  0  0] ; 
line_dir=[0  0  1] ; 
bvec=latVec*bvec  n ' ; 


%unit  burgers  vector  /  slip  direction 
%disocation  line  direction 
%burgers  vector 


%  RDX  properties  (MPa)  Table  3.6,  p79  of  Munday  thesis  (UMd  2011)  and 
%  Munday  et  al .  (JPCB  2010) 

%  Voigt  Stiffness  Coefficients  from  Munday  Dissertation 

Cl l=25e3 ; 

C22=23 . 8e3 ; 

C33=23 . 4e3 ; 

C4  4=3 . Ie3 ; 

C55=5 . 2e3; 

C6  6=7 . 7e3 ; 

C23=8 . 8e3 ; 

C13=7 . 6e3; 

C12=10 . 6e3 ; 
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■END  MANUAL  INPUT 


%  Voight  notation 


Cv= [Cll 

C12 

C13 

0 

0 

0;  . 

C12 

C22 

C23 

0 

0 

0;  . 

C13 

C23 

C33 

0 

0 

0;  . 

0 

0 

0 

C44 

0 

0;  . 

0 

0 

0 

0 

C55 

0;  . 

0 

0 

0 

0 

0 

C6  6]  ; 

S=inv (Cv) ; 

%  Find  average  isotropic  material  Properties 

%  Voigt  Average  -  uniform  strain,  over-estimate  of  stresses 

%  bulk  modulus 

Bv=0  ; 

for  i=l : 3 


(Get  from  C) 


Bv=Bv+sum (Cv ( i , 1 : 3 ) ) ; 
end 

Bv=Bv/ 9 ; 

%  shear  modulus 

Gv= ( Cv (1,1) +Cv (2,2) +Cv (3,3) ) /15- 

(Cv (1,2) +Cv (2,3) +Cv (1,3) ) /15+ (Cv ( 4 , 4 ) +Cv ( 5 , 5 ) +Cv ( 6 , 6 ) ) /5; 

%  Youngs  Modulus  &  Poissons  Ratio 
Ev= ( 9*Bv*Gv) / ( 3*Bv+Gv) ; 
nuv= (3*Bv-2*Gv) / 2 / (3*Bv+Gv) ; 

%  Reuss  Average  -  uniform  stres,  under-estimate  of  stressses  (Get  from  S) 

Br=0; 
for  i=l : 3 


Br=Br+sum ( S ( i, 1 : 3) ) ; 
end 

Br=l/Br; 

%  Shear  Modulus 

Gr=4/15*  (S  (1, 1)  +S  (2, 2)  +S  (3,3)  )  -4/15*  (S  (1,2)  +S  (2, 3)  +S  (1,3)  )  +1/5*  (S  (4,  4)  +S  (5, 5)  +S  (6,  6)  ) 
Gr=l/Gr; 

%  Youngs  Modulus  &  Poissons  Ratio 

Er= ( 9*Br *Gr ) / (3*Br+Gr) ; 
nur= (3*Br-2*Gr) / 2/ (3*Br+Gr) ; 


fprintf ( ' \n\n ISOTROPIC  PROPERTIES : \n ' ) 
fprintf ( '  Voigt  Avg  Properties:  Ev=%6.2f, 
\n ' , Ev, nuv, Gv, Bv) 

fprintf ( '  Reuss  Avg  Properties:  Er=%6.2f, 
\n ' , Er, nur, Gr, Br ) 


nuv=%6 . 2f , 
nur=%6 . 2f , 


Gv=%6.2f  , 
Gr=%6.2f  , 


Bv=%6 . 2f 
Br=%6 . 2f 


%  get  fourth  order  C 
C=zeros (3, 3, 3, 3)  ; 
for  i=l : 3 
for  j=l:3 
for  k=l : 3 
for  1=1:3 
m=0; 
n=0  ; 
if  i==j 
m=i  ; 

elseif  (  (i==2  &&  j==3)||(i==3 
m=4  ; 

elseif  ( (i==l  &&  j==3) | [ (i==3 
m=5  ; 

elseif  (  (i==l  &&  j==2)||(i==2 
m=6  ; 
end 

if  k==l 
n=l  ; 

elseif  (  (k==2  &&  1==3)  |  |  (k==3 


&&  j  ==2 ) ) 
&&  j  ==1 ) ) 
&  &  j  ==1 ) ) 


&&  1=2 )  ) 
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n=4  ; 


elseif 
n=5  ; 

\ — 1 

II 

II 

2^ 

&&  1==3 ) | | 

(k==3 

&&  1==1) ) 

elseif 

( ( k==l 

CM 

II 

II 

i — 1 

<-3 

ua 

( k==2 

\ — 1 

II 

II 

i — 1 

<-3 

ua 

n=6; 

end 

C  (i,  j  ,  k,  1)  =Cv  (m,  n)  ; 
end 
end 
end 
end 

% - CALCULATE  KK - 

%permutation  or  Levi-Cevita  3x3x3  tensor  (using  linear  indexing  of  matrix) 
lc  =  zeros  (3,3,3) ; 
lc ( [8  12  22] )  =  1; 
lc ( [6  16  20] )  =  -1; 

%  Integral  for  equation  2.2 

t=line_dir /norm ( line_dir ) ;  %disocation  line  direction  --  eq  2.3 
r=norm (t) ; 

theta=atan2 ( t ( 2 ) , t  ( 1 ) ) ; 
phi=acos (t  (3 ) /r) ; 

a= [sin (theta)  -cos (theta)  0];  %eq  2.5 

d= [cos  (phi) *cos (theta)  cos  (phi) *sin (theta)  -sin (phi)];  %  eq  2.5 

%Integral  Range:  psi  ranges  0  to  pi 

psi_l=0; 

psi_2=pi ; 

Ninc=200;  %  #  of  midpoint  summation  intervals 

psi_inc= (psi_2-psi_l ) /Nine; 
intzM=zeros (3, 3, 3, 3) ; 
intzMiso=zeros  (3, 3, 3, 3) ; 

LAM3=zeros (3,3) ; 

for  inc=l:Ninc  %midpoint  integration  loop  for  intgral  of  christophel  matrix 

psi=psi_l+ (inc-1 ) *psi_inc; 

z(l:3)  =a*cos (psi) +d*sin (psi) ;  %eq  2.4 

dz (1 : 3) =-a*sin (psi ) +d*cos (psi) ;  %eq  2.4 

%  Christoffel  stiffness  --  C_ijrs  z_r  z_s 
M=zeros  (3,3) ; 
for  i=l : 3 
for  r=l : 3 
for  j  =1 : 3 
for  s=l:3 

M ( i , r ) =M ( i , r ) +C ( i , j , r , s ) * z  ( j ) * z  ( s ) ; 
end 
end 
end 
end 

%midpoint  integration  of  Christophel  terms  for  eq  2.2 

invM=inv (M) ; 
for  s=l:3 
for  n=l : 3 
for  i=l:3 
for  r=l:3 

intzM(s,n,i,r)=intzM(s,n,i,r)  +  ( z ( s ) *dz (n) *invM ( i , r ) ) *psi_inc ; 
end 
end 
end 
end 
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end 


%  SOLVE  FOR  ANISOTROPIC  DISLOCATION  COMPLIANCE  OR  ORIENTATION  PREFACTOR  TENSORS 
%  this  is  the  inverse  of  LAMBDA  given  by  Sun  &  Beltz  (JMPS  1994)  etc. 
KK_BA=zeros (3, 3) ;  %  Barnett  &  Asaro  JMPS  v.20  1972 

for  m=l : 3 
for  g=l:3 
for  p=l : 3 
for  j  =1 : 3 
for  w=l : 3 
for  n=l : 3 
for  i=l : 3 
for  r=l : 3 
for  s=l : 3 

%  --  Equation  2 . 2 

KK_BA (m, g) =KK_BA (m, g) +1/ 8/piA2*lc  (p, j , w) *t  ( j )  .  .  . 

* (C (n,g, i,p) *C (w,m, r, s) +C (n,m, i,p) *C (w, g, r,  s) )  . . . 
*intzM ( s, n, i , r ) ; 

end 

end 

end 

end 

end 

end 

end 

end 

end 


%%%% - OUTPUT  RESULTS - 


%  Dislocation  energy 

Eaniso=bvec ' * (KK_BA) *bvec  *  le6  *  le9;  %  eq  2.1 
cos_theta=dot (bvec, line_dir/norm (line_dir) ) /norm (bvec) ; 

Eiso  v=l/4 /pi*Gv* ( l-nuv*cos  thetaA2 )/( 1-nuv) * (norm (bvec) ) A2* le6* le9 ; %eq  2.6 

Eiso_r=l/4 /pi*Gr * ( l-nur*cos_thetaA2 ) / ( 1-nur ) * (norm (bvec) ) A2*le6*le9; 

fprintf ( ' \n\nENERGY  PRELOG  FACTOR  ( JmA-l*le-9) :  \n\n') 

fprintf ( '  Anisotropic  =  %-6g\n ' , Eaniso) 

fprintf ( '  Isotropic  Voigt  =  %-6g\n ' , Eiso_v) 

fprintf ( '  Isotropic  Reuss  =  %-6g\n ' , Eiso_r) 


%  splitting  distance  between  two  edge  partials 

d_aniso= (4*pi*bvec ' * (KK_BA) *bvec)  /2/pi/ (y_sf*le-6) *lelO;  %eq  3.17 

d_isov=  (Gv/ (1-nuv) ) /2/pi/ (y_sf*le-6) * (norm(bvec) ) A2*lel0; 

d_isor=  (Gv/ (1-nur) ) /2/pi/ (y_sf*le-6) * (norm(bvec) ) A2*lel0; 

fprintf (' \n\nEDGE  PARTIAL  SPLITTING  DISTANCE  (A):  \n\n ' ) 

fprintf ( '  Anisotropic  =  %- 6g\n ' , d_aniso) 

fprintf ( '  Isotropic  Voigt  =  %- 6g\n ' , d_isov) 

fprintf ( '  Isotropic  Reuss  =  %- 6g\n ' , d_isor ) 

"6 

%  NUCLEATION  FOR  KII  LOADING  FOR  EMISSION  OF  EDGE  DISLOCATION 
%  equation  3 . 8 

fprintf (' \n\n  KII  LOAD  FACOTRS  FOR  FIRST  PARTIALS  (MPa  sqrt (m) ) :  \n\n') 
fprintf ( '  Isotropic  Voigt  =  %-6g\n sqrt (2*Gv/ (1-nuv) *y_us*le-6) )  %eq 
fprintf ( '  Isotropic  Reuss  =  %-6g\n sqrt (2*Gr/ (1-nur ) *y_us*le-6) ) 
fprintf ( '  Anisotropic  =  %- 6g\n ' , sqrt ( 8 *pi*KK_BA ( 1 , 1 ) *y_us* le-6 ) ) 

%  equation  3.21 

fprintf (' \n\n  KII  LOAD  FACOTRS  FOR  SECOND  PARTIAL:  \n\n') 
fprintf ( '  Isotropic  Voigt  =  %-6g\n sqrt (2*Gv/ (1-nuv) * (y_us-y_sf) *le-6) ) 
fprintf ( '  Isotropic  Reuss  =  %-6g\n sqrt (2*Gr/ (1-nur )* (y_us-y_sf) *le-6) ) 
fprintf ( '  Anisotropic  =  %- 6g\n ' , sqrt ( 8 *pi*KK_BA ( 1 , 1 ) * ( y_us-y_sf ) *le- 6 ) ) 
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